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ABSTRACT 

Diffusive acceleration at collisionless shock waves remains one of the most 
promising acceleration mechanisms for the description of the origin of cosmic 
rays at all energies. A crucial ingredient to be taken into account is the reaction 
of accelerated particles on the shock, which in turn determines the efficiency of 
the process. Here we propose a semi-analytical kinetic method that allows us to 
calculate the shock modification induced by accelerated particles together with 
the efficiency for particle acceleration and the spectra of accelerated particles. 
The shock modification is calculated for arbitrary environment parameters 
(Mach number, maximum momentum, density) and for arbitrary diffusion 
properties of the medium. Several dependences of the diffusion coefficient on 
particle momentum and location are considered to assess the goodness of the 
method. 
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1 INTRODUCTION 



Most scenarios for the origin of cosmic rays rely upon the acceleration of charged particles 
in the presence of shock waves, developed in sources such as supernova remnants, active 
galaxies, planetary shocks, gamma ray bursts and many others. The ba sic feature s of th e 
acceleration process have b e en highligh ted in the pioneering papers bv iKrvmskiil ()1977l ): 
Blandford fc Ostrikerl (jl978[ ): lBelll (jl978T ) in the context of the so-called test particle assump- 
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tion. Several excellent reviews ([Drurvl (|1983f ) ; iBlandford fc Eichlerl (jl987h ; 
( 19911) : lMalkov fc Drurvl 1)20011 )) discuss in detail the many problems that are still open con- 
cerning the acceleration process. Among these, a fundamental one is the limited applicability 
of the results found within the test particles approach. In most scenarios for the origin of 
either galactic or extra-galactic cosmic rays, in fact, an appreciable fraction of the kinetic 
energy crossing the shock needs to be transfered to accelerated particles. This need contra- 
dicts the very assumption that the accelerated particles are test particles, unable to exert 
any dynamical reaction onto the shocked fluid. The well known result that the spectrum 
of the accelerated particles is a power law with slope nearly independent of the detailed 
properties of the system (e.g. diffusion coefficient) holds only within the context of this test 
particle approximation. Relaxing this assumption leads to the modification of the shock 
by the accelerated particles, a pheno menon that ha s received m uch attention in the con- 



text of the so-called two-fluid mode 



(12221); 



Malkov. Diamond & Volkl (200 



s (IDrurv Volkl (I1980i Il98lh l. kinetic models (jMalkov 



d: IniasJ J2002I kooi 



and numerica l approaches , 



both Monte Carlo and other s imula ti on procedures (Jones & 



Ellison Mo 



'1997 



2005); 



jius fc Paschmannl (ll990T) 



Kane. 



ones 



lllison. Baring fc Jones! ()1995 . 



1996): 



llison ill 991): iBelll (Il 987): 



Kane fc Jones 



fc Gieselerl (j2002J)). For an accurate recent review see the work 



bv iMalkov fc Drurv (2001), from which the reader can see the weak and strong points of 
each approach. The present paper illustrates a kinetic analytical approach, which provides 
the exact solution for the spectrum of accelerated particles and shock modification in a very 
general situation in which the diffusion properties of the medium are arbitrary. The problem 
is reduced to solving an integral-differential equation, which easily leads to the required so- 
lution. For th e injection of particles at t he sh ock surface we implement the recipe previously 



presented by 



Blasi. Gabici &: Vannonil (J2005). In all the cases that we considered we never 



find evidence for multiple solutions. The method we propose is of general validity, in that 
it can be used for an arbitrary momentum dependence of the diffusion coefficient and for 
diffusion properties (related to the magnetization properties of the medium) that can change 
in an arbitrary way with the spatial location in the fluid. 



2 CALCULATIONS 

The equations for the conservation of the mass and momentum fluxes between upstream 
infinity and a point x in the upstream region can be written as: 
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Po«o = p(x)u(x), (1) 

Po ul + P gfl = p{x)u{xf + P g (x) + P C r{x) , (2) 

where p, u and P g are the gas density, velocity and pressure (the corresponding quantities 
at upstream infinity are indicated with the index 0). The pressure of accelerated particles is 
defined as 



Pcr{x) 



Pmax 



dp 4:iip 3 v(p)f(x,p), 



(3) 



Pinj 



and f(x,p) is the distribution function of accelerated particles. Here pi n j and p m ax are the 
injection and maximum momentum. The function / vanishes at upstream infinity, which 
implies that there are no cosmic rays infinitely distant from the shock in the upstream 
region. The distribution function satisfies the following transport equation in the reference 
frame of the shock: 



d_ 

dx 



d 

D(x,p)—f(x,p) 



u df(x,p) + } L fdu s 
dx 3 \ dx . 



p^ + Q(*, p ) = o. 



(4) 



The x axis is oriented from upstream infinity (x = — oo) to downstream infinity (x = +oo), 
with the shock located at x — 0. The injection is introduced here through the function 
Q(x,p). The diffusion properties are described by the arbitrary function D(x,p), depending 
on both momentum and space. In previous approaches restrictive assumptions on the dif- 
fusion coefficient were always adopted in order to facilitate the path to the solution. The 

solution / can be written in the following implicit form: 

u(x') 



f(x,p) = exp 



Up) + 



l r° dx' 



3 Jx D(x',p) 



exp 



dx" 



U[X 



11 \ 



D(x",p) 



dx' 



p2 Qp 



D(x',p) 



dx' 



du 
dx" 



f(x" lP )p 3 



In the case of a spatially constant diffusion coefficient, as shown bv iMalkovl (1993 ), a very 



(5) 



good approximation to the solution is found in the form f(x, p) = fo(p) exp — J° dx'u(x') 



with q(p) = — d d , and fo(p) = f(x = 0,p) the distribution function at the shock. We found 



dlnp 

that the similar form 



f(z,p) = /o(p)exp 



q(p) f° , u(x') 
dx 



D(x\p) 



(6) 



represents a very good approximation for the case of diffusion coefficients with arbitrary 
spatial dependence (see Sec. 3). We adopt therefore this functional form in our calculations, 
although it is not strictly required, in the sense that we could well use the complete solution, 
Eq.0 
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It was shown by 

3 Rtot 



Blasil ((20021 ) that the function fo(p) can be written in general as 



/o(p) 



7]n 



R tot U(p)-lJ iirpl 



exp 



p dp 7 ZRtotUW) 
m] p' Rt ot U (p')-l 



Here we introduced the function Z7(p) = u p /u , with 



Wr 



Ml - 



/o(p) 



dx(du/dx)f(x,p) 



(7) 



(8) 



where «i is the fluid velocity immediately upstream (at x = ). We used Q(x,p) = 
vn ? as ^ Ux 8{p — p in j)5(x), with n gas i = rioRtot/Rsub the gas density immediately upstream 
(x = 0~) and rj the fraction of the particles crossing the shock which are going to take part 
in the acceleration process. In Eq. Q we also introduced the two quantities R su b = Ui/u 2 
(compression factor at the subshock) and R to t = Uq/u2 (total compression factor). The two 
compression factors, assuming, for simp l icity, that the heating is only adiabatic, are related 
through the following expression ((Blasil ((2002)): 



R 



tot 



Mr 



( 7s + l)R] 3 ub - ( 7s - 1)R 



7 9 +l 
sub 



7ff + l 



(9) 



where M is the Mach number of the fluid at upstream infinity and 7 5 is the ratio of 
specific heats for the fluid. The parameter t] in Eq. [7| contains the very important in- 
formation about the injection of particles from the thermal bath. Following the work of 

(2005(), we relate f] to the compression factor at the subshock as: 

(10) 



Blasi. Gabici fc Vannoni 



V 



37T 1 / 



m {R sub - l)£ 3 e-« 



Here £ is a parameter that identifies the injection momentum as a multiple of the momentum 
of the thermal particles in the downstream section (p ini = Hvthp)- This recipe is i nspired to 



the thermal leak age model originally presented b y 



previous work by 



Ellison. Jones fc Eichler ( 1981 ): Ellison! ( 198lL 



Giesele r. Jones fc Kand (1200011 ( see a. 



Ellison fc Eichlerl (jl984h l 



so 



The parameter £ is supposed to contain the information about the microscopic structure 
of the shock. For collisionless shock waves the thickness of the shock is expected to be of 
the order of the Larmor radius of the thermal particles in the shock vicinity, which is not 
a very well defined concept because of the violent fluctuations in the electromagnetic fields 
within the shock. A simp l e argu ment can be used to infer that £ is of the order of 2 — 4 
( Blasi. Gabici fc Vannonil ((20051 ) ). For the numerical calculations that follow we always use 
£ = 3.5, that allows for only a fraction of order 10 -4 of the particles crossing the shock to 
be injected in the accelerator. 
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Eq. El for the conservation of the momentum flux, once normalized to PqUq, is easily 



transformed to 



£ c (x) = 1 + 



U{x) 



1 



;U(x) 



lg Ml K ' l9 Ml 
where £ c (x) = Pcr{x)/poUq and U(x) 
(Eq. [HJ, we can also write: 



-7 9 



(11) 

u(x)/uq. In terms of the distribution function 



Air 
3p Mo 



P rn a x 



dp p 3 v(p)f (p)exp 



x p (x',p) 



(12) 



where for simplicity we introduced x p (x,p) = ■ 
By differentiating Eq. ^] with respect to x we obtain 



^ = -X(x)Ux)U(x), 
ax 



where 



X{x) =< 1 / x p >e c 



(x',p) 



(13) 



(14) 



g:; x dp P 3 v(p)fo(p) exp [- £ dx'^^ )} 
and U (x) is expressed as a function of £ c (x) through Eq. HU 

Finally, after integration by parts of Eq. |H1 one is able to express U(p) in terms of an 
integration involving U (x) alone: 



U(p) 



dx U (x) 



Xp (x, p) 



exp 



- , V u< - x "> 



(15) 



Xp {x , p) ^ 

which allows to easily calculate fo(p) through Eq. UJ 

Eqs. [H] and E3 can be solved by iteration in the following way: for a fixed value of 
the compression factor at the subshock, R su b, the value of the dimensionless velocity at 
the shock is calculated as U(0) = R su b/Rtot- The corresponding pressure in the form of 
accelerated particles is given by Eq. [TT]as £ c (0) = 1 + ^jjj - - (-fpf) 79 • This is 
used as a boundary condition for Eq. El where the functions U(x) and X(x) (and therefore 
fo(p)) on the right hand side at the k th step of iteration are taken as the functions at the 
step (k — 1). In this way the solution of Eq. ^2 at the step k is simply 



er J w = e c (o)ex P 



dx'X^ix'jU^ix') 



(16) 



with the correct limits when x — > and x — > — oo. At each step of iteration the functions 
U(x), fo(p), X(x) are recalculated (through Eq. HH Eqs. Hoi and 171 and Eq. HH respectively), 
until convergence is reached. The solution of this set of equations, however, is also a solution 
of our physical problem only if the pressure in the form of accelerated particles as given by 
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p/ mc 



Figure 1. Upper panel: Spectra of accelerated particles at the location of the shock for Mo = 4 (dotted line), 10 (short-dashed 
line), 50 (dash-dotted line), 100 (dash-dot-dot-dotted line), 300 (long-dashed line) and 500 (solid line). Lower panel: momentum 
dependent slope for the same values of Mach numbers. In both panels we used £ = 3.5 and p m ax = 10 5 me. 

Eq. ^2 coincides with that calculated by using the final fo(p) in Eq. ED This occurs only for 
one specific value of R su b, which fully determines the solution of our problem. 

3 RESULTS 

The computational method illustrated in the previous section is very fast and allows one 
to determine the solution (namely the velocity and density profiles in the precursor, the 
density of accelerated particles as a function of momentum and location in the upstream 
fluid and all the thermodynamic quantities related to the gas) for an arbitrary choice of 
the diffusion coefficient and for any values of the environmental parameters (Mach number, 
density, maximum momentum). 

In Fig. [U we illustrate the spectra (upper panel) and slopes (lower panel) as a function 
of momentum for the following values of the Mach number: Mq = 4 (dotted line), 10 (short- 
dashed line), 50 (dash-dotted line), 100 (dash-dot-dot-dotted line), 300 (long-dashed line) 
and 500 (solid line). The distribution functions are multiplied by p A to emphasize the concave 
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Mach Number Mo 




Rtot 


&(0) 


Pin J 




V 




4 


3.19 


3.57 


0.1 


0.035 


3.4 


x 10" 


-4 


10 


3.413 


6.57 


0.47 


0.02 


3.7 


x 10" 


-4 


50 


3.27 


23.18 


0.85 


0.005 


3.5 


x 10" 


-4 


100 


3.21 


39.76 


0.91 


0.0032 


3.4 


x 10" 


-4 


300 


3.19 


91.06 


0.96 


0.0014 


3.4 


x 10" 


-4 


500 


3.29 


129.57 


0.97 


0.001 


3.5 


x 10" 


-4 



Table 1. Shock modification for different Mach numbers. 

shape of the modified spectra. All the curves refer to p m ax = 10 5 in units of mc. The most 
evident aspect of shock modification, found in all previous calculations, is here confirmed: 
the shock modification is enhanced when the Mach number of the shock increases. The 
spectrum is flatter at high momenta as confirmed by the lower panel of Fig. Q and easily 
understood in terms of the large values of the total compression factor (see Table EJ). 

For strongly modified shocks, the slo pe becomes even flatter tha n p~ 3 - 5 at high momenta, 



as also found in numerical simulations (jBerezhko fc Ellison! (jl999l ) and refs. therein) 1 . In 



these conditions, most energy is channeled in the highest energy part of the spectrum. At 
lower energies on the other hand, the spectrum is steeper than that predicted by linear 
theory, as a natural consequence of the lower compression at the subshock for strongly mod- 
ified shocks. For the parameters adopted here, the energy saturation (namely £ c (0) ~ 1) is 
achieved for Mach numbers around 100, as demonstrated by the fact that the corresponding 
curves in the upper panel of Fig. ^ have roughly the same height (namely the same energy 
content). On the other hand, different modifications result in different compressions at the 
subshock and therefore different injection momenta. This is illustrated in Fig. and Table 01 
In particular in Table 01 we list the values of the compressione factors, dimensionless cosmic 
ray pressure at the shock, injection momentum and fraction of accelerated particles for the 
same values of M used to obtain the curves in Fig. [TJ 

In Fig. El we illustrate the results of our method for different choices of the momentum de- 
pendence of the diffusion coefficient. We consider three cases: 1) Bohm diffusion, Ds{p) oc p; 
2) Kraichnan diffusion, Dx r (p) ocp 1 ^ 2 ; 3) Kolmogorov diffusion, D^oiip) ocp 1 ^ 3 (relativistic 
scalings). For illustrative purposes, we choose to calculate the spectrum of accelerated par- 
ticles and the shock modification for Mo = 100 and p m ax = 10 5 mc. The resulting spectrum 
is shown in the upper panel of Fig. |2J for Bohm diffusion (solid line), Kraichnan diffusion 
(dotted line) and Kolmogorov diffusion (dashed line). The general tendency is that the sat- 



1 We remind that in other semi-analytical calculations (e.g. iMalkovl I1997I) ') the asymptotic spectrum for Pi„j <p < p m ax has 
slope 3.5 
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Figure 2. Upper panel: Spectra of accelerated particles at the location of the shock for Mo = 100, p m ax = 10 5 mc and for a 
Bohm diffusion coefficient (solid line), Kraichnan diffusion coefficient (dotted line) and Kolmogorov diffusion coefficient (dashed 
line). Lower panel: Distribution of the pressure in the form of accelerated particles, normalized to the ram pressure (£, c {x) as 
defined by Eqs. 1111121 . for the same three cases. The spatial coordinate is in units of x* = Ds(p m ax)/«o, with Db the Bohm 
diffusion coefficient. 



uration phenomenon occurs at somewhat lower Mach numbers for diffusion coefficients that 
depend more weakly on momentum. The lower panel in Fig. El illustrates the spatial distri- 
bution of energy in the accelerated particles (£ c (x)), where the spatial coordinate is chosen 
in such a way that x — 1 in the point = DB(Pmax)/uo- Clearly the particles with the 
maximum momentum diffuse on shorter spatial scales than for diffusion coefficients with 
weaker momentum dependence. 

The power of the computational method in being suitable for treating arbitrary depen- 
dences of the diffusion coefficient on momentum and spatial coordinates is further demon- 
strated in Fig. 01 where we show how the solutions change when the diffusion coefficient is 
allowed to vary in space. For illustrative purposes we consider the case of a Bohm diffusion 
coefficient with D B (p,x) oc p (constant in space) and D B (p,x) oc p/p(x), where p(x) is the 
gas density at the position x, self-consistently calculated by using the conservation laws. 
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Figure 3. Upper panel: Spectra of the accelerated particles for spatially constant Bohm diffusion (thin lines) and for Bohm 
diffusion with D(p,x) oc p/p(x). The different line-types refer, for each case, to x = (solid line), x = 10~ 7 x» (dot-dashed 
line), x = 10 -4 x» (short-dashed line) and x = 0.1 x* (long-dashed line). The dotted lines neighbouring each curve refer to the 
distribution functions computed by using in the right-hand-side of Eq. |F]the solution obtained with our method (see text for 
details). Lower panel: S; c (x) and U(x) (solid and dashed lines, respectively) for the case of spatially constant Bohm diffusion 
(thin lines) and for D(p,x) oc p/p(x) (thick lines). The spatial coordinate is again in units of x* defined as for Fig. 121 



The latter dependence is representative of the case of a magnetic field frozen in the plasma 
flowing in the upstream section. 

In the upper panel of Fig.Elwe plot the spectrum of the accelerated particles for spatially 
constant Bohm diffusion (thin curves) and for Ds{p,x) oc p/p(x) (thick lines). The different 
line-types refer to spectra at the different spatial locations: x = (solid line), x = 10~ 7 x* 
(dot-dashed line), x = lO -4 ^* (short-dashed line) and x = 0.1a;* (long-dashed line), where 
x* is defined as above, i.e. a;* = DB(p m ax)/uo with Db referring to the spatially constant 
Bohm diffusion coefficient. In the lower panel we plot £ c (a;) (solid lines) and U(x) (dashed 
lines) for the same two cases, identified by the different thickness of the lines. 

In order to assess the goodness of our approximate solution (Eq. EJ) we computed the 
right-hand-side of Eq. El by using the functions U(x) and f(x,p) found with our method. 
The correction is found to be non- negligible only in the exponentially decreasing parts of the 
spectrum (see dotted lines in Fig. EJ), which contains negligible energy and hardly leads to 
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any observable features. On this basis, we conclude that Eq. El is an excellent approximation 
to the solution for diffusion coefficients with arbitrary spatial and momentum dependences. 
The solutions obtaine d with this method are remarkably similar to those obtained with 



approximate methods by 



Blasi (2002 



2004), for the case of Bohm diffusion with no spa- 



tial dependence. The discrepancies with such previous treatments are expected and indeed 
appear for increasingly weaker dependences of the diffusion coefficient on the particles' mo- 
mentum, and in general when a spatial dependence of the diffusion properties is added (these 
aspects will be discussed in an up coming det ai led p aper) . The results also compare w ell with 



a previous method proposed bv lMalkovl (jl997l ) and lMalkov. Diamond fc Vol k (2000), where 



the contribution of gas pressure was neglected and no recipe for injection was adopted. 

The most important property of the method here described, however, is the fact that 
it appears to be the first that allows to take into account the spatial dependence of the 
diffusion coefficient. The importance of being able to deal with arbitrary diffusion prop- 
erties is highlighted by the following considerations. First, particle acceleration at shocks 
i s expected to be efficient o nly if the turbulenc e responsible for diffusion is self-generated 



(jLagage fc Cesarskvl (|l983h : lLucek fc Belli (l2000M Belll j2QQ4j )), and in this case the diffusion 
coefficient is necessarily dependent upon both momentum and space in a complex manner. 
Moreover, the appearance of a maximum momentum is indeed due to the fact that at some 
distance from the shock diffusion becomes uneffective and particles are no longer trapped in 
the shock vicinity. Since the shock modification depends in a crucial way on the value of the 
maximum momentum, it is clear that a careful calculation of the shock modification should 
be able to account for these phenomena. 
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